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1.   INTRODUCTION 

Recent  interest  in  the  numerical  integration  of  stiff 
ordinary  differential  equations  has  led  to  a  search  for  methods  of 
higher  order  than  those  currently  known.   This  paper  describes  a 
computer  graphic  technique  which  was  developed  to  make  a  systematic 
search  possible.   New  high  order  methods,  found  using  this  technique, 
will  be  presented. 

Stiff  equations  are  characterized  by  the  presence  of  widely 
varying  time  constants.   Terms  corresponding  to  very  small  time 
constants  decay  rapidly  and  are  of  no  interest  because  of  their  small 
magnitude.   Conventional  methods  are  unstable  if  the  step  size  used 
is  much  greater  than  the  smallest  time  constant.   Stiffly  stable 
methods  allow  a  much  larger  step  size  while  maintaining  stability. 

Known  stiffly  stable  methods  are  predictor-corrector 
methods.   The  corrector  may  be  written 

(1"1)         \J»tViJ«-i,'"noVi 

where  a  ^  0,  |a  |  +  -g  J  >  0.   The  discretization  error  of  the  method 

K.  U        (J 

is  defined  as  the  difference  between  the  value  of  y  calculated  from 

n 

(l-l)  and  the  value  of  the  exact  solution  y(t  ).   Define 

e  -  y  -  y(t  )  =  y  -  F 

n   Jn  J      n    ^n    n 


Then  the  errors  e  obey  the  difference  equation 


(1-2) 


I      (a.  +  hXB.)  e  _.  +  T(F,  t  ,  h)  =  0 
.  -   1      i   n-k+i      *  n' 
i=0 


where  T  is  the  truncation  error  for  the  exact  solution.   Convergence 

requirements  make  it  necessary  to  hound  the  solution  of  (1-2).   This 

[2] 

requires    the  stability  of  the  homogeneous  difference  equation 

obtained  from  (1-2)  with  T  set  to  zero.   This  equation  is  stable  if 

and  only  if  the  modulus  of  no  root  of  the  polynomial  equation 

(1-3)  k 

£   (a.  +  hAB. )  r   =  0 
i-0  1 

exceeds  1,  and  roots  of  modulus  1  are  simple. 

The  order  p  of  a  multistep  method  may  be  defined  as  the 

highest  degree  of  polynomial  functions  y  for  which  the  truncation 

[3] 
error  T(y,  t  ,  h)  vanishes  identically  in  t  and  h.     The  truncation 

error  is  then  directly  proportional  to  h   ,  so  p  should  be  as  large 

as  possible.   Dahlquist    has  shown  that  if  (l-3)  is  stable  for  all 

A  with  negative  real  parts,  then  the  maximum  order  of  the  method  is  two, 


Gear,  C.  W. ,  "The  Numerical  Integration  of  Stiff  Differential 
Equations,"  University  of  Illinois,  Department  of  Computer 
Science,  Report  No.  221,  p.  k. 

[2] 

Henrici,  P.,  Discrete  Variable  Methods  in  Ordinary  Differential 

Equations,  John  Wiley  and  Sons,  Inc.,  New  York,  1962,  pp.  218-220. 

[3] 

Gear,  C.  W.  ,  ojo_.  cit.  ,  p.  5. 

[h] 

Dahlquist,  G. ,  "A  Special  Stability  Problem  for  Linear  Multistep 
Methods,"  BIT  3,  1963,  pp.  27-^3. 


Gear    has  relaxed  these  requirements  and  defined  those  necessary 
for  stiff  stability.   The  requirements  are  shown  in  Figure  1. 

hA  Plane 
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-0 


Figure  1, 


Since  accuracy  is  usually  defined  in  terms  of  order,  where  order  is 
defined  as  above,  we  again  see  that  a  large  p  is  desirable.   Previous 
results  show  the  existence  of  stiffly  stable  methods  of  order  as  high 
as  six  for  proper  D,  0,  and  a.   The  object  of  the  current  investigation 
was  to  find  new  methods  of  order  greater  than  six. 


[5] 


Gear,  C.  W.  ,  op_.  cit .  ,  pp.  5-7. 


2.   METHOD  OF  TESTING 

k  k 

Define  p(?)  =  E  a.^1  and  o(c)  =  £   8.51.   Then  (1-3) 
i=0  X  i=0  x 

may  be  rewritten 

(2-1)  pU)  +  hXa(e)  =  0 

Letting  hX  ■*■  °°,  the  roots  of  (2-1)  approach  those  of 

o(£)  =  0.   This  implies  that  the  roots  of  o(£)  must  also  lie  inside 

the  unit  circle  or  be  on  the  unit  circle  and  simple,  and  that  a(c)  is 

of  degree  k.   Given  any  o(£),  the  unique  method  of  degree  k  may  be 

\  f>  1 
found  as  follows.   For  a  method  of  degree  k,  we  have 

(2_2)  P(0  +  log(0  0(e)  ~  c(C-Dk+1  as^l 

k  , 

Set  £  -  1  =  x.   From  (2-2),  p(l+x)  =  -  [  Z   6.  (l+x)J]  log(l+x), 

j=0  J 
truncated  to  terms  in  x  and  lower. 

rn 

Let  ■$■  =    [aQ,  a±9 ,   afc]  ,   then 


L    ■*      Henrici,   P.,   op_.    cit.  ,   pp.    226-228. 
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The  last  matrix  consists  of  the  coefficients  of  the  logarithm  series 

for  (l+x).   Each  matrix  in  the  inner  sum  of  matrices  corresponds  to  a 

k 
multiplication  by  (l+x)   in   £   8.  (l+x)  .   The  first  matrix  provides 

j=0  J 


the  conversion  back  to  £  =  l+x. 


We  want  hX  values  such  that  (2-l)  has  roots  inside  the 
unit  circle  or  on  the  unit  circle  and  simple.   This  region  is  bounded 

by  the  locus  of  -  ^y^f  in  the  hX  plane  for  £  =  e10 ,  0e[0,2tt].   If 

a(c)  is  stable,  the  method  will  be  stable  at  hX  =  infinity,  so  by 

continuity  arguments  any  region  connected  to  hX  =  infinity  will  be 

[7] 

stable.     In  order  to  obtain  values  for  plotting  this  locus,  we 

let  y(Q)  =  -  (fy  witn  £  =  el"  and  calculate  y(0)  for  0  =  0(^)tt, 

where  M  is  chosen  to  provide  a  small  enough  increment  to  obtain  a 
suitable  number  of  points.   The  calculation  of  y(Q)    is  as  follows 

y(e)  =  -  »*& 


o(C) 

kiG 
k 

10 

0     kiG 

3ke         +... 

...+  B/9  +  B0 

Letting   R  =  £  a   cos(jQ)   ,  S  =   E   a  sin(jO) 

j=0  J  j=l  J 

k  k 

A  =  T.  B.  cos(j0)   ,  B  =  I   B,  sin(jO) 


we  get 


j=0  J  j=l  j 


fn\  R  +  iS 

Y(0)  =  "  aTTb 

(RA  +  SB)  +  i(SA  -  RB) 

2    2 
A  +  B 


w   1  +       (RA  +  SB)    ,   _    (SA  -  RB) 
Now  let  x  =  -  -^ — ^—  and  y  =  -  — 5 p— 

A  +  B  A  +  B 


[7] 

L,J   Gear,  C.  W.  ,  op_.  cit.  ,  pp.  9-10. 


The  desired  locus  is  obtained  by  plotting  the  curve 
passing  through  the  x,y  values  calculated  for  ©  varying  over  the  range 
0  to  f,  giving  the  lower  half,  and  then  rotating  this  curve  180*  to 

obtain  the  upper  half,  since  the  locus  is  symmetric  about  the  real 

v 
axis.   Figure  2  shows  a  sample  plot  for  o(C)  =  K     with  k  =  6. 


-«.o 


Figure  2. 


Given  a  k  and  the  coefficients  of  the  polynomial  o(£), 
the  stiff  stability  of  the  unique  method  of  degree  k  may  thus  he 
investigated  by  visually  examining  the  plotted  locus.   The 
computer  graphic  approach  used  in  the  current  study  takes  advantage 
of  this  fact.   A  plot  of  the  locus  is  created  on  a  CRT  display  and 
an  interactive  system  is  provided  for  manipulating  k  and  the 
coefficients  of  a(£).   The  great  advantage  is  that  of  the  speed 
of  obtaining  a  new  plot.   Conventional  plotting  techniques  are 
too  slow  to  make  a  study  of  a  large  number  of  methods  possible  in 
a  reasonable  length  of  time. 

While  speed  was  a  major  goal  in  designing  the  graphic 
system,  it  should  be  noted  that  implementation  of  the  system  is 
on  a  DEC  PDP-8  computer  with  an  attached  Model  338  CRT  display. 
The  PDP-8  is  a  12  bit  per  word,  fixed  point  binary  machine  whose 
basic  arithmetic  capability  is  two's  complement  fixed  point  addition. 
The  necessary  floating  point  arithmetic  is  thus  accomplished  entirely 
by  means  of  software.   It  was  therefore  desirable  to  minimize  the 
number  of  calculations  necessary  to  obtain  the  points  used  in 
plotting  the  locus . 

The  calculation  of  p-  given  in  (2-3)  may  be  optimized  in 
terms  of  fewest  floating  point  operations  as  follows.   Given  k 
and  6  ,  B-,  » ,  B,  generate  the  Pascal  triangle  matrix  A,  where 

A. .  =  (.),  or  as  shown  in  (2— U )  and  where  the  element  in  the  upper 

left  corner  of  A  is  indexed  A   and  so  on. 


(2-U) 


A  = 


1111 

12   3 

1  3 

1 


(k) 


k 
1 


To  generate  the  matrix  B  representing  Z   8.  (l+x)  ,  for 

j  =  0,  l*,..t  k  multiply  the  j-th  column  of  A  by  the  corresponding  6. 

J 


1 


r  "\ 

i 


6o  {      }h< 


0 


0 
v  J 


r   ^ 

i 


)'"  ^    ( 


(k) 
\2j 


(k) 

k 

1 


10 


Sum  across  the  rows  to  get  the  column  matrix  M, 
(2-5) 


M  = 


l      6. 
j=0  J 


k 

£  J  6. 
3-1    J 


M 


0 


M 


"k 


Since  A  is  an  upper  triangular  matrix,  indexing  on  the  implemented 
versions  of  the  above  operations  is  constructed  in  such  a  manner 
that  additions  and  multiplications  are  performed  only  on  upper 
triangular  elements.   In  terms  of  (2-5),  B  is  now  given  below  by 


B  = 


M, 


Ml  Mo 

M2  Mx 


M. 


\     \-l     \-2    •  •  •  Ml  M0 


11 


The  large  savings  in  number  of  multiplications  and 
additions  over  straightforward  calculation  as  in  (2-3)  is  obvious. 
In  addition,  a  significant  savings  in  storage  during  calculation 
and  for  B  itself  is  realized.   Because  of  the  nature  of  B,  only 
M  need  be  retained.   A  transformation  on  the  indexing  function 
on  B  may  be  used  to  obtain  corresponding  elements  from  M. 

The  first  matrix  in  (2-3)  is  obtained  by  negating  each 
element  of  A  the  sum  of  whose  subscripts  is  odd.   Call  the  new 
matrix  A' .   The  next  step  in  the  calculation  of  p-  is  the  matrix 
multiplication 

A*  x  B  =  C 
Taking  advantage  of  the  fact  that  A1  is  upper  triangular  and  B 
lower  triangular,  and  using  M  to  obtain  the  desired  elements  of 
B,  the  most  efficient  method  of  obtaining  C  is  summarized  in  the 
following  section  of  ALGOL  program 

for  i  :=  0  step  1  until  k  do 

for  ]  ;=  0  step  1  until  k  do 
begin  sum  :=  0.0; 

for  n  :=  i  step  1  until  k  do 
if  j  <  n  then 

sum  :=  sum  +  A(i ,n)*M(n-j ) ; 
C(i,j)  :=  sum; 
end; 
Multiplication  of  the  matrix  C  by  the  last  matrix  in  (2-3)  then 
gives  us  -  jy  =  -  [oQ,  o  , ,  akJ  . 


12 


In  order  to  increase  the  efficiency  of  the  calculation 

of  y(Q),  tables  of  values  of  sin(O)  and  cos(O)  for  ©  =  0(— )n  are 

stored  permanently.   A  table  look-up  technique  is  employed  to 
obtain  the  values  of  sin(jO)  and  cos(jO)  as  needed.   The  technique 
works  in  the  following  manner.   Each  0  will  be  of  the  form 

ItTT 

— ,  k  =  0,  1,  2, ,  M.   The  table  entries  for  sin  and  cos 

may  be  thought  of  as  being  numbered  from  0  to  M.   Thus  the  value 

of  cos(— )  is  found  in  entry  number  k  in  the  cos  table.   To  obtain 

values  of  sin(jO)  and  cos(jO)  we  note  that  j0  =  *•— : — — .   If  we 

form  i  =  ( j •k)mod(M) ,  then  the  desired  value  will  be  *  entry 
number  i  in  the  sin  or  cos  tables  respectively.   The  correct 
sign  may  be  determined  by  counting  the  number  of  times  M  must 
be  subtracted  from  j*k  to  get  i.   An  odd  number  indicates  that 

a  sign  change  is  necessary.   This  is  because  jQ  =  *•— : - —  falls 

in  a  quadrant  diagonally  opposite  to  that  —  falls  in.   As  an 

example  let  M  =  100,  k  =  30,  and  j  =  h,    so  that  j0  =      =  — — . 

Then  i  =  20  and  one  subtraction  of  100  is  necessary  to  get  i,  so 
the  desired  sin  or  cos  is  minus  entry  20  in  the  respective  tables. 
It  should  be  obvious  from  the  above  that  incrementing 

over  0  =  0(— V  simply  corresponds  to  incrementing  k  over  k  =  0(l)M. 
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Values  of  sin(jO)  and  cos(j0)  for  use  in  calculating  any  particular 
point  on  the  locus  are  thus  obtained  using  an  i  based  on  j  and  the 
current  value  of  k  to  index  the  tables. 


Ik 


3.   SYSTEM  DESCRIPTION 

The  graphic  system  was  designed  to  make  a  rapid  search 
for  new  stiffly  stable  methods  possible.   Easy  to  use  interactive 
capabilities  were  thus  included  to  control  the  display  and  the  basic 
input  parameters  k  and  o(£).   Implementation  is  on  a  PDP-8  with  four 
U096  word  memory  banks.   The  attached  Model  338  display  runs  off 
display  files  constructed  in  the  computer  memory.   An  on-line 
teletype  and  a  bank  of  12  push-buttons  are  used  for  interactive 
communication  by  the  current  system. 

Computer  storage  allocation  is  shown  in  Figure  3. 


Bank :  0 
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DISPLAY 
FILE  B 

FLOATING 

POINT 

SOFTWARE 
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COORDINATES 
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;  ....  _ 

Figure  3. 


The  majority  of  the  system  programs  in  bank  0  are  concerned  with 
the  calculation  of  -p  and  y{&)   using  the  methods  described  in  the 
previous  section.   The  others  calculate  an  appropriate  display  file 
for  each  locus. 
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The  display  is  based  on  a  user  defined  display  window 
specified  in  terms  of  left  (XT )  and  right  (X^)  boundaries  relative 

to  the  real  (X)  axis,  and  upper  (Y.J  and  lower  (Y  )  boundaries 

relative  to  the  imaginary  (Y)  axis.   Only  the  section  of  the  locus 
the  coordinates  of  which  fall  within  these  boundaries  is  displayed. 
Axes  are  displayed  as  appropriate  to  the  specified  window.   Each 
display  is  scaled  to  cover  the  entire  usable  screen  area.   The  338 
screen  size  is  102U  units  by  102U  units,  where  a  unit  is  an 
intensifiable  point,  and  points  are  numbered  0  through  1023.   A 
coordinate  x  will  be  displayed  if  X  <  x  <  X  ,  and  the  calculation 

x  -  X 

*102U 


performs  the  required  scaling.   A  similar  calculation  produces  the 
desired  y'  for  YT  '-   y  <  Y  .   The  display  consists  of  vectors 

connecting  consecutive  points  (x' ,y' )  on  the  locus.   If  the  0 
increment  used  in  calculating  y(Q)»  and  thus  the  (x' ,y' )  coordinates, 
is  small  enough,  the  locus  will  appear  as  a  smooth  curve. 

In  order  to  increase  the  efficiency  of  the  display  file 
building,  an  appropriate  switch  is  set  whenever  a  point  on  the  locus 
first  falls  above,  below,  left,  or  right  of  the  display  window.   When 
each  new  point  thereafter  is  tested,  the  switch  on  causes  the 
corresponding  out  of  window  condition  to  be  tested  immediately.   Thus 
only  a  single  test,  rather  than  a  four-way  test,  is  necessary  on  each 
new  point  as  long  as  the  condition  occurs.  When  the  locus  re-enters 
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the  window,  a  non-intensified  vector  is  drawn  from  the  last  point 
in  the  window  to  the  new  point,  and  the  switch  is  turned  off. 

The  most  obvious  requirement  of  the  system  is  to  allow  a 
user  to  input  a  k,  the  coefficients  8  ,  B,  , ,  6,  of  the 

polynomial  a(£),  and  the  display  window  boundaries,  from  which  a 

display  of  the  -  / ~\    locus  will  be  created.   Since  it  is  much  too 
a(fcj 

cumbersome  to  respecify  all  of  this  data  each  time  a  specific 

parameter  is  to  be  changed,  interactive  facilities  for  selectively 

changing  the  parameters  have  been  included  in  this  system.   The 

amount  of  recalculation  necessary  to  create  a  new  display  based  on 

the  change  varies  significantly  depending  on  the  nature  of  the  change. 

Two  display  files  are  maintained  so  that  the  old  display  is  left  on 

until  the  new  display  file,  resulting  from  the  latest  parameter 

change,  has  been  built.   Display  control  is  then  transferred  to  the 

new  file. 

User  intentions  with  respect  to  specific  changes  are 
indicated  by  the  use  of  the  push-buttons,  each  of  which  has  an  implied 
meaning.   Messages  asking 'for  new  parameters  are  printed  on  the  teletype 
as  required,  after  which  the  user  inputs  the  appropriate  data.   Push- 
buttons are  associated  with  the  following  actions. 

1.   Respecify  k,  SQ,  6^ ,  6k  and  the  3^,  X^  Yy,  YL 

display  window  boundaries.   A  new  locus  is 
automatically  displayed. 
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2.  Change  g. .   User  specifies  i,  then  new  3..   No  new 

locus  is  calculated.   This  allows  the  user  to  change 
more  than  one  3.  before  getting  a  new  display. 

3.  Calculate  a  new  locus  and  create  a  new  display  based 
on  one  or  more  3.  changes. 

h.      Specify  an  increment  for  8. .  User  specifies  i  and 

the  increment  amount  by  which  8.  is  to  be  incremented. 

1 

No  incrementing  is  done  here. 

5.  Increment  8.  based  on  the  data  specified  in  the  last 

use  of  the  button  associated  with  action  h.      A  new 
locus  is  calculated  and  displayed  using  the  new  B.. 

This  is  separated  from  action  k   so  that  a  3.  may  be 

incremented  without  requiring  new  teletype  input  each 
time  it  is  incremented. 

6.  Respecify  the  display  window  boundaries.   A  new  display 
is  created  based  on  the  new  window. 

7.  Double  display  window  size  and  create  new  display.   Since 
no  teletype  input  is  required,  this  enables  a  very  fast 
change  in  the  display  when  a  larger  viewing  area 

is  desired. 

8.  Halve  display  window  size  and  create  new  display. 

9.  List  the  (x,y)  coordinates  of  the  points  on  the 
currently  displayed  locus. 
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10.  List  the  a    ,   a    , a     coefficients   of  the 

U    _L  K 

polynomial  p(£)  calculated  in  connection  with  the 
currently  displayed  locus . 

11.  List  the  current  6^,  3n  , ,  6,  coefficients.   This 

0   1         k 

is  particularly  useful  after  incrementing  or  changing 

several  3. . 

l 
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k.      TESTS  AND  RESULTS 

Using  the  system  described  in  the  previous  section,  a 
systematic  search  for  stiffly  stable  methods  of  order  greater  than 
six  has  been  conducted.   Suitable  coefficients  for  a(£)  have  been 
found  which  indicate  the  existence  of  new  methods  of  order  as  high 
as  eight. 

In  addition  to  the  stability  requirement,  the  relationship 
between  the  3's  and  the  number  of  items  which  must  be  saved  to 
provide  information  about  the  history  of  the  solution  influenced 
the  choice  of  polynomials  o(£)  tested.   The  corrector  formula  (l-l) 
can  be  written  more  generally  as 

a  y  (m+l)  =  -  E  a    y      -  h  E1  6    y'   (m+l) 
k  yn+l  k-j  yn-j+l   n  .=Q  \-j   yn-j+l 

By  defining  v  as 

^=    [V   W '  yn-k+l'   h  yn>   h  K-V '   h  Ck^l1 

r  8 1 

and  defining  suitable  A,   &_,   and  F,   Gear         has   shown  that  the  numerical 
predictor-corrector  process  may  be  written  as 

(m+l)  (m)        .      /      (m)  x 

yn+l         =  yn+l     +^F(^n+l    >   Xn+1'  h) 


[8] 


Gear,  C.  W. ,  "The  Numerical  Solution  of  Ordinary  Differential 
Equations  of  Various  Orders,"  Argonne  National  Laboratory  Report 
ANL  7126,  Argonne,  Illinois,  January  1966,  pp.  7-9. 


20 

For  general  o(£)  as  defined  in  this  study,  we  would  have  k  =  k  in 

[k-l).      However,  if  3.  =  0,  i  =  0,  1, ,  j ,  then  k  =  k-j-1. 

Since  y  contains  information  to  be  saved  at  step  n,  it  is  obvious 

from  (k-2)   that  small  k  require  the  saving  of  less  information  at 

each  step  of  the  corrector  iteration. 

In  order  for  (k-2)   to  make  sense,  we  must  have  k  >_  1. 

The  least  amount  of  solution  history  data  must  be  saved  when  k  =  1. 

Looking  at  (4-l),  this  corresponds  to  allowing  only  B  and  B    to 

have  non-zero  values.   Initial  search  efforts  were  thus  directed 
to  methods  involving  a(£)  satisfying  this  condition. 

Since  methods  of  order  up  to  six  were  known,  the  most 
immediate  goal  of  the  current  investigation  was  to  find  a  stiffly 
stable  method  of  order  seven.   For  convenience,  B  was  set  equal 

to  one  throughout  the  study,  since  the  coefficients  of  o(£)  may 
always  be  scaled  in  such  a  manner  as  to  make  this  hold.  Fixing 
B  to  one,  and  imposing  the  above  restriction,  made  it  possible  to 

study  the  most  optimal  possibilities  of  order  seven  by  varying  only  B/-. 

Stability  requires  the  roots  of  o(C)  to  lie  within  the  unit  circle,  or 
be  onthe  unit  circle  and  simple.   Another  restriction  on  a(£)  is  that 
it  have  no  common  factor  with  p(£).   Since  p(^)  will  always  have  a 
root  at  +1,  a(C)  must  not  have  a  root  at  +1.   In  terms  of  the  above 
restriction,  with  B  =1,  this  implies  that  B    must  not  equal  -1.   For 

k  K—-L 

optimal  methods  of  order  seven,  possible  B/-  are  thus  limited  to  the 
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range  -1  <  3,-  f_  1.   By  setting  3,-  equal  to  -.99,  a  locus  was 

obtained  which  indicated  the  existence  of  suitable  parameters 
D  and  0  for  a  stiffly  stable  method.   Thus  the  unique  method 
determined  by  the  polynomial  a(£)  =  £  -  •  99£  was  found  to  be  a 
stiffly  stable  method  of  order  seven.   As  3/-  was  incremented  toward 

+1,  the  parameter  0  was  seen  to  decrease  toward  zero  until  3,-  passed 

-.8,  at  which  time  a  suitable  0  ceased  to  exist.   Details  of  this 
and  other  new  methods  are  summarized  in  the  Appendix. 

Using  this  same  technique,  a  similar  study  of  the  most 
optimal  possibilities  for  orders  up  to  15  was  made.   No  locus  was 
found  which  indicated  the  existence  of  a  stiffly  stable  method. 
Results  indicate  that  no  parameter  0  exists  for  methods  of  these 
orders  determined  by  stable  a(£)  satisfying  the  3, 


k-2   Mk-3 


&     =  0  condition. 


By  removing  the  restriction  on  non-zero  3's,  less  optimal 
possibilities  for  methods  of  order  greater  than  six  were  investigated. 
Efforts  were  concentrated  on  methods  of  order  seven  and  eight.   For 
these  orders,  each  coefficient  3.   i  =  0,  1, ,  k-1  was  incremented 

over  the  range  -1  <  3.  £  1  with  3,  =  1  and  3.  =  0,  j^i,j^k. 

1  k         J 

For  k  =  7,  loci  indicating  stiffly  stable  methods  were  found  for  3> 


in  the  range  .6  <_  3>  <.  .8.   Thus,  for  example,  the  unique 


method 


7     h 
determined  by  the  polynomial  o(£)  =  £  +  .75   is  stiffly  stable.   No 

stable  methods  were  found  when  this  approach  was  applied  to  methods 
of  order  eight. 
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5 .   CONCLUSIONS 

Although  no  simple  method  for  theoretically  determining 
stiffly  stable  numerical  integration  methods  is  known,  the  computer 
graphic  system  described  here  has  enabled  the  efficient  examining  of 
large  numbers  of  possible  methods.   Systematic  search  procedures 
have  led  to  the  discovery  of  stiffly  stable  methods  of  orders  seven 
and  eight. 

Defining  optimality  in  terms  of  least  amount  of  solution 
history  data  to  be  saved  at  each  step  of  the  corrector  iteration, 
the  optimal  possibilities  for  orders  seven  through  fifteen  were 
tested,  with  negative  results  for  all  orders  except  seven.   Results 
of  the  tests  indicate  that  no  optimal  methods  of  order  greater  than 
seven  exist. 

Removing  the  optimality  restriction  allows  the  unique 
method  determined  by  any  stable  a(c)  to  be  a  valid  candidate  for 
satisfying  stiff  stability  requirements.   No  general  results  can  be 
given,  since  the  number  of  potential  methods  excludes  the  possibility 
of  any  type  of  exhaustive  testing.   Testing  of  large  classes  of  methods 
tends  to  indicate  the  scarcity  of  stiffly  stable  methods  of  orders 
greater  than  six.   Methods  of  orders  seven  and  eight  were  found,  however, 
while  examining  some  of  the  more  likely  classes  of  possible  methods.   The 
greatest  advantages  of  the  graphic  system  seem  to  be  the  ability  to 
rapidly  accept  or  reject  whole  classes  of  a(£),  and  the  ability  to  narrow 
in  on  a  method  in  a  promising  region  by  interactively  manipulating  the 

coefficients  of  a(c)  while  watching  the  changes  in  the  locus  of  -  /    \ . 
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This   section  -will  present   details   of  the  new  methods 
of  order  greater  than  six.      Two  methods   of  order  seven  and  one  of 
order  eight  will  be   given. 

The   coefficients   for  an  optimal  method  of  order  seven  are 

iO  1  2  3  k  5  6 


3. 

i 

0 

0 

0 

0 

0                 0 

-.99 

1 

a. 

l 

0.166 

-1.361+ 

U.9^2 

-10.1+00 

lU.lUl  -13.U70 

8.1+35 

-2.1+51 

The  locus  of  -  / A   is  given  in  Figure  A-l.   From  Figure  A-l  we  see 

that  -6.1  is  an  approximate  upper  hound  on  the  parameter  D  and  that 
.5  is  an  approximate  upper  hound  on  the  parameter  0. 
Another  method  of  order  seven  is  given  by 


0 


.7 


a.        0.138       -1.120         3.990       -8.050       11.1+92  -10.920  7.070     -2.560 


The  locus   of  -    t't  for  this  method  is   given  in  Figure  A-2.      An 
approximate  upper  bound  for  D  is   -12.1+  and  for  0,   .25. 
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A  method  of  order  eight  is  specified  by 


i  0      1 

2 

3       h 

5 

6 

7 

8 

S.  0      0 

i 

0 

0       0 

0 

.81 

-1.8 

1 

a.    -.±62      1.U892 
1 

-6.129 

1U.890  -23.763 

26.587 

-21.070 

IO.636 

-2.U78 

The  locus  of  -  ~JTT  ^s  given  i-n  Figure  A-3.   The  approximate  upper 
bound  for  D  is  -6.9  and  for  0,  .25. 


27 


|    I     I    I    I     |     I     I     I     I     |     I    I 


100 


Figure  A-l 


28 


- 10.0/     =  ; 


Figure  A-2. 
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Rigure  A- 3. 
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